Dissipative dynamics of a qubit coupled to a nonlinear oscillator 
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We consider the dissipative dynamics of a qubit coupled to a nonlinear oscillator (NO) embedded 
in an Ohmic environment. By treating the nonlinearity up to first order and applying Van Vleck 
perturbation theory up to second order in the qubit-NO coupling, we derive an analytical expression 
for the eigenstates and eigenfunctions of the coupled qubit-NO system beyond the rotating wave 
approximation. In the regime of weak coupling to the thermal bath, analytical expressions for 
the time evolution of the qubit's populations are derived: they describe a multiplicity of damped 
oscillations superposed to a complex relaxation part toward thermal equilibrium. The long-time 
dynamics is characterized by a single relaxation rate, which is maximal when the qubit is tuned to 
one of the resonances with the nonlinear oscillator. 

PACS numbers: 03.67.Lx,03.65.Yz,05.45.-a,85.25.-j 



I. INTRODUCTION 



Coupling a two-level system (TLS) to a harmonic os- 
cillator has attracted a lot of attention in various fields 
of physics. Examples are two-level quantum dots in pho- 
tonic crystal nanocavities 1 '- 2 , a quantum dot exciton in a 
microcavity^, or single atoms with a large dipole moment 
interacting with photons in a microwave cavitjA Within 
the framework of quantum computation two prominent 
solid-state realizations of a qubit-oscillator system are 
found: a Cooper-pair boy£&iZ& coupled to a transmis- 
sion line resonato r 9 ' 10 ! 11 ' 12 ! 13 ' 14 and the Josephson flux 
qubil^ read-out by a DC-SQUIDikiL 1 ^. The Cooper- 
pair box setup has been used to perform non-demolition 
measurements or to transfer information between qubits 
via the transmission line resonato r 9 ' 12 ' 19 ' 20 ' 21 ' 22 . In the 
second experimental realization the flux qubit is usually 
read-out via a damped DC-SQUID, which acts as a lin- 
ear or nonlinear oscillator. A non-demolition read-out 
procedure, based on the measurement of the Josephson 
inductance, is given by Lupa§cu et ali^ 3 .. 
At present the effort to exploit the nonlinearity of a qubit 
read-out device, for example, a DC-SQUID or a Joseph- 
son bifurcation amplifier (JBA) 24 ' 25 , is growing, as non- 
linear effects lead to advantages in various measurement 
schemes and to new physical observations. For exam- 
ple, the qubit read-out can be optimized by using the 
SQUID in the nonlinear regime as a bifurcation amplifier 
leading to fast read-out with high fidelit y 26 ' 27 . Second, 
the bifurcation allows for a higher sensitivity when deter- 
mining the qubit states and, due to the nonlinear Joseph- 
son inductance, a high quality factor for the resonance is 
achieved^. However, the nonlinear regime also provides 
new channels of relaxation^ 7 -. Moreover there are recent 
experiments embedding a micromechanical resonator in a 
nonlinear DC-SQUID, which is strongly damped to avoid 
bistability, to acquire cooling and squeezing of the res- 
onator modes and to achieve quantum-limited position 
detection^. Such a composed system can then also be 
coupled to a qubit. Besides these examples a SQUID 
which is embedded into a cavity— can be used as a bi- 



furcation amplifier in its nonlinear regime. 
All these approaches rely in principle on treating the 
SQUID as a classical nonlinear system. To our knowl- 
edge there has been to date no experimental realization 
of a SQUID in the nonlinear quantum regime. 
From the theoretical point of view nonlinear quantum os- 
cillators have been predominantly studied within the con- 
text of the quantum Duffing oscillator model^i^Li 3 ^ 3 -! 3 ^, 
where the oscillator is subject to an external ac driving 
force. Strikingly, the response of the Duffing oscillator 
displays antiresonant dips and resonant peaks depending 
on the frequency of the driving field^ .. The antireso- 
nances persist in the presence of a weak Ohmic bath; 
for high damping the nonlinear response of the oscilla- 
tor resembles the one of a linear oscillator at a shifted 
frequenc y 31 ' 32 ^ 3 -. 

Despite the numerous theoretical works on coupled qubit- 
linear oscillator systems^^^^i 3 ^^^^^ the case of 
a TLS- Josephson bifurcation amplifier system has been 
addressed only very recently by Nakano et alJ*. Here 
we study the SQUID as a nonlinear, undriven oscilla- 
tor acting as a read-out device for a qubit. We consider 
weak nonlinearities such that the corresponding linear 
system can be retained at any step of our calculation. 
With the help of Van Vleck perturbation theory in the 
TLS-oscillator coupling g we determine the eigenstates 
and spectrum of the coupled system and the correspond- 
ing dynamics in analytic form. Thus we can quantita- 
tively characterize the influence of the coupling g and of 
the nonlinearity on the dynamics of the composed sys- 
tem. The overall effects of the nonlinearity are the fol- 
lowing: (i) a shift of the transition frequencies to higher 
values compared to the linear case; (ii) the amplitudes 
associated to the transition frequencies are modified. In 
particular the vacuum Rabi splitting is decreased by the 
interplay of coupling and nonlinearity. To account for 
dissipative effects we add a weak Ohmic environment. 
Then the dynamics of the reduced density matrix of the 
composed system can be described in terms of a set of 
coupled differential equations for its matrix elements in 
the energy basis (Bloch-Redfield equations). We discuss 
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a partial secular approximation (PSA) to those equations 
as well as two more stringent approximations, the full 
secular approximation in the low temperature approxi- 
mation (LTA) and the smallest eigenvalue approximation 
(SEA) accounting for the long time dynamics. All these 
three approximation schemes allow for analytical solu- 
tion of the dynamics of the TLS, which we compare with 
predictions obtained by numerically solving the Bloch- 
Redfield equations. It turns out that the most accurate 
PSA should be used when investigating strong nonlin- 
earities. The long-time approximation enables us never- 
theless to extract the correct relaxation rate within the 
regime of validity of our perturbative approach. The pa- 
per is organized as follows: In section [ill we introduce 
the model with the relevant dynamical quantities. In 
section ITTIl the energy spectrum and the dynamics of the 
non-dissipative coupled system is investigated. Section 
IIVI addresses the dissipative effects, while in section [V] 
results are represented. In section [VI] conclusions are 
drawn. 



II. THE MODEL 

A. Qubit-nonlinear oscillator-bath system 

In this section we consider a TLS coupled to a nonlin- 
ear oscillator, which itself is coupled to an Ohmic bath. 
This model mimics, e.g., the situation of a flux qubit, 
made of three Josephson junctions, which is coupled in- 
ductively to a damped DC-SQUIL^I. The qubit with 
its two logical states, the clockwise and counterclock- 
wise currents, represents a two-level system. Because the 
SQUID itself is coupled to an environment, it transfers 
environmental influences which lead to the dissipation in 
the qubit. Hence the total Hamiltonian reads: 



T~C — Wtls-no + ^no-b + 7~Cb, 



(1) 



with Wtls-no describing the coupled TLS-nonlinear os- 
cillator system, while Hno-b and Hb are the coupling 
between the oscillator and bath and the bath Hamilto- 
nian, respectively. For later convenience we write 



^tls-no — Htls + ~M.NO +Wln 



(2) 



with coupling Hamiltonian Wmt- 



1, Two-level system 



First we consider the Hamiltonian of the TLS, 



(3) 



represented in the localized basis {\L),\R}}^, corre- 
sponding to clockwise and counterclockwise currents, re- 
spectively, in the superconducting ring. The 



are the corresponding Pauli matrices. The energy bias e 
can be tuned for a superconducting flux qubit by applica- 
tion of an external flux <j> oxt and vanishes at the so-called 
degeneracy points. 

For e 3> Ao, where Ao is the tunneling amplitude, the 
states \L) and \R) are eigenstates of the TLS, while at 
the degeneracy point the eigenstates |g), |e) are given 
by symmetric and antisymmetric superpositions, respec- 
tively, of the two logical states. In general the states \R) 
and \L) become in the energy basis: 



\R) =cos(e/2)|g)+ S m(e/2)|e), 
\L) =-sin(e/2)|g) + cos(e/2)|e> 



(4) 



with tan 9 = — A /e and — 5 < 9 < |. Moreover in 

this basis the TLS Hamiltonian is: Htls — — ^2^ z ^ 
where a z is the Pauli matrix in the energy basis and 
hAb = Ti\J E 1 + Aq is the energy splitting. 



2. Nonlinear oscillator 

The Hamiltonian for the nonlinear oscillator is com- 
posed of a linear harmonic oscillator modified with a 
quartic term in the position operator, 



Hno = hn] + -(B + B^)\ 



(5) 



where j = B^B is the occupation number operator of 
the linear oscillator and B and B^ are the corresponding 
annihilation and creation operators. In the following we 
restrict to the case of hard nonlincarities, i.e., a > 0. 
Using time-independent perturbation theory we consider 
small nonlinearities a <C htt and evaluate the eigenvalues 
£j and eigenfunctions \j) of ([5]) to lowest order in the 
nonlinearity, 



Sj := mj+-aj(j + l), j 



o,. 



, oo 



(6) 



li) := |j)o + a°]|j-2)o + 4 j) b- + 2)o+ (7) 
a^|j-4) + 4 3) |j+4) , 

where |)o denotes the eigenstate of the corresponding lin- 
ear oscillator. The expansion coefficients for the jth state 
of the nonlinear oscillator are given by: 



,0') 



(j) 



,0) 



v /(j-3)(j-2)(j-1)jq 
V(j + l)(j + 2)(j+3)(j + 4)a 

2m 



(8) 



(j + §) y/(j + l)(j + 2)a 



2nn 
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Error 


q/M = 10~ 3 


a/Ml = 0.01 


a/Ml = 0.02 


£r«(l) 


3 ■ ltr 3 


0.03 


0.06 


Er {2 \l) 


2.06 ■ 10~ 5 


2.06 ■ 10" 3 


8.25 • 10" 3 


Er {1 \2) 


4.5 • 10 -3 


0.045 


0.09 


Er [2) {2) 


3.84 ■ 10" 5 


3.8 • 10" 3 


0.015 


Er w (3) 


6 ■ 10 -3 


0.06 


0.12 


Er< a >(3) 


6.56 • 1(T 5 


6.56- 10" 3 


0.026 


£r (1) (4) 


7.5 • 10" 3 


0.075 


0.15 


£r< 2 '(4) 


1.02 ■ 10~ 4 


1.02 ■ 10 -2 


0.041 


£r (1) (5) 


9 ■ 10" 3 


9 ■ 10~ 2 


0.18 


£r (2) (5) 


1.46 • 10" 4 


1.46 • 10" 2 


0.058 



TABLE I: Error estimation for different values of the nonlin- 
earity for the six lowest levels. 



We notice that two arbitrary eigenstates \j), \k) are or- 
thonormal up to first order in the nonlinearity. 
Perturbation theory for a nonlinear oscillator has to be 
elaborated carefully. Due to the special form of the non- 
linear term, proportional to (B + B^) 4 the energy correc- 
tions acquire a strong level dependence: Sj = §cr/(j+l) 
for the first, see Eq. ©, and sf> = g^a 2 (-34j 3 - 
51j 2 — 59j — 21) for the second order. Depending on the 
actual level number the second order can be as large as 
the first order for fixed nonlinearity. To avoid this, one 
has to choose the nonlinearity parameter a such that the 
oscillator levels under consideration are well represented 
by the first order result. The error done by disregarding 
the nth order perturbation theory is estimated in the fol- 
lowing by introducing Er^ n \j) — \£^\/£j for different 
nonlinearities (see Table [I]). Taking only first order per- 
turbation theory into account, the error is determined by 
Er^ (j max ), where j ma x is the highest level under consid- 
eration. The error made by using first order perturbation 
theory is in case of a/Ml — 0.02 around 6% for the j = 5 
level. 

Finally we consider a coupling Hamiltonian of the form: 



Hi„t = hga z (B + B^) 



(9) 



This kind of coupling arises due to the inductive coupling 
of the TLS to the SQUID^. 



3. Harmonic bath 

Following Caldeira and Leggetl4£, we model the envi- 
ronmental influences originating from the circuitry sur- 
rounding the qubit and the oscillator as a bath of har- 
monic oscillators being coupled bilinearly to the nonlin- 
ear oscillator. Thus, the environment is described by 
TLb = Y]h Mjkb\pk and the interaction Hamiltonian is 



Hno-b = {B^ + B)J2Hv k {bl + b k ) + (B^+BfJ2H 



(10) 



The operators b\ and b k are the creation and annihilation 
operators, respectively, for the fcth bath oscillator, ujk is 
its frequency, and v% gives the coupling strength. The 
whole bath can be described by its spectral density, which 
we consider to be Ohmic, 

GohmM = ^ ~ u k) = KW, (11) 



where k is a dimensionless coupling strength. 



B. Population difference 

We wish to describe the dynamics P(t) of the TLS 
described by the population difference 



P(t) = Tr TLS {cr z p re d(0} 

= (R|p red (i)|R)-(L|p red (t)|L) 



(12) 



between the \R) and \L) states of the qubit. The reduced 
density matrix of the TLS, 



p xed (t) = Tr NO Tr B {IF(i)} = Tt NO {p(t)}, 



(13) 



is found after tracing out the oscillator and bath de- 
grees of freedom from the total density matrix W(t) = 
exp~Ti nt W(0) exp"R wt . For vanishing nonlinearities it is 
possible to map the problem described by the Hamilto- 
nian in equation (fTJ) onto a spin-boson models with an 
effective peaked spectral density depending on the cou- 
pling (7, the frequency f2, and the damping strength k. 
This mapping hence allows the evaluation of the popula- 
tion difference P{t) of the TLS using standard approxi- 
mations developed for the spin-boson mode l 39 ' 40 ^. This 
mapping, however, does no longer hold true in the non- 
linear oscillator case. Hence in this work we consider the 
TLS and the nonlinear oscillator as central quantum sys- 
tem and describe dissipative effects by solving the Bloch- 
Redfield master equations for the reduced density matrix 
p(t) = TrB{W(t)} of the qubit-NO system. In a second 
step we perform the trace over the NO degrees of freedom 
to obtain the reduced dynamics of the TLS. An expres- 
sion for P(t) is then given in terms of diagonal and off- 
diagonal elements of pit) in the 7Ytls-no Hamiltonian's 
eigenbasis It reads^: 



where 



Pnn(t) = 



^2 jcose 



> 2 - (je\n) 2 



(14) 



(15) 



2 sin 9 ( jg | n) (je\n)} p nn (t), 



2 |cos9 



(jg|n)(m|jg) - (je\n)(m\je) 



sin0 



(je\n}(m\jg) + (je\m)(n\jg) 
Re{p nm (t)}, 
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and p nm (t) — (n\p(t)\m) . The TLS-NO eigenstates are 
derived in the next section. 



III. ENERGY SPECTRUM AND DYNAMICS 
OF THE NON-DISSIPATIVE TLS-NO SYSTEM 

In the following we derive the eigenenergies and 
eigenstates of the unperturbed TLS-NO Hamiltonian 
Htls-no using Van Vleck perturbation theor y 50 ' 51 . 
This approach allows us to deal with spectra containing 
almost exactly degenerate levels organized in manifolds 
(here doublets), as it is the case if the TLS and nonlinear 
oscillator are close to resonance, ~ fi, and the cou- 
pling g is small compared to the energy separation of the 
manifolds. 



A. Energy spectrum 

The eigenstates of the uncoupled TLS-NO system 
Hamiltonian H are {\j) <g) |g); \j) ® |e)} = {|jg); \je)}. 
The associated energies are depicted by the dotted lines 
in figure [1] At the resonance condition of the TLS with 
two neighboring nonlinear oscillator levels, 

HQ = hA b -3a(j + 1), (16) 

where j denotes the lower oscillator level involved, the 
states |(j + l)g) and |je) are exactly degenerate except 
for the ground state |0g). For finite coupling the full 
Hamiltonian Htls-no acquires in the basis {|ig);|je)} 
the form 

Wtls-no = Ho + Hint (17) 

= — + m 3 + + ^ + 

^ (ea z - AqCTz) (B + 5 f ) . 

To find the eigenvalues of the Hamiltonian 7Ytls-nOj w ^ 
treat 7ii n t cx g as a small perturbation, which is satisfied 
for g <C Ah, f2. Using Van Vleck perturbation theor y 50 ' 51 
we can construct an effective Hamiltonian by applying 
an unitary transformation to Wtls-nOj 

Weff = exp(iS , )7i: TL s_ N oexp(-iS'). (18) 

Ties has the same eigenvalues as Htls-no but does not 
involve matrix elements connecting states which are far 
away from degeneracy. Consequently it is block-diagonal 
with all quasi-degenerate energy levels being in one com- 
mon block. Because the quasi-degenerate states form 
doublets, each block of H e fi is given by a 2 x 2 matrix. 
The latter can be diagonalized easily. To calculate S and 
7i e ff we write both as a power series up to first order in 
the nonlinearity a and up to second order in the coupling 
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FIG. 1: Energy spectrum of the coupled qubit-nonlinear- 
oscillator system versus the linear oscillator frequency £1 (in 
units of the TLS tunneling splitting Ao). Solid lines show 
the energy levels for the five lowest energy states (|0), |1), 
|2), |3>, |4>) with the TLS-NO coupling being switched on, 
g = O.I8A0, and for finite nonlinearity, a = 0.02ftAo. The 
TLS is unbiased, e — 0. The energy levels for the uncoupled 
case are given by the dotted lines. Due to the non-equidistant 
level spacing of the nonlinear oscillator the resonance condi- 
tion (crossing of dotted lines), given in equation (|16[) . is differ- 
ent for each doublet. This causes a shift of the exact crossings 
with respect to the linear case at zero coupling to lower fre- 
quencies. For finite coupling the spectrum exhibits avoided 
crossings around resonance, whereas it approaches the uncou- 
pled case away from resonance. 



9, 

S = S^+S^+S^+0(a 2 ,g% (19) 
H eS = W<°> + W$ + H<$ + 0(a 2 , g 3 ), (20) 

where exp(iS^) = 1. The upper index in the above 
equation denotes the actual order in g. Consequently 
in the following we assume that a/HQ, ~ g 2 /£l 2 < 1. 

To calculate S^ 1 / 2 ) and 7Y^ /2) we use both that Tl Q g acts 
only inside a manifold and that S has no matrix elements 
within a manifold. The general formulas are found e.g. 

m 49.50.51 

The results for the effective Hamiltonian and the trans- 
formation matrix are given in the appendix [X] The non- 
vanishing matrix elements of the effective Hamiltonian, 
apart from the zeroth-order contributions in g, are 

(n cS f\ , =-m^n 1 ( j ) = HA( j ), (21) 

^ /ie;(j+l)g A 6 

and 

(h cS ) {2) = fi[WiO\n)-w<, (?',«)], (22) 
(h cS ) (2) = h\Wi{j,a) + w a (j + i,n)]. (23) 
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We used as abbreviation 



ni(j) = 1 + 



1 - 



V7TT 
3q 



2M] 



O' + i) 



V7TT j 

+ 0(a 2 ), (24) 



Wo(j.") 



A^(A b + f!) 



(26) 

Therefore the effective Hamiltonian acquires in first order 
in the nonlinearity and in second order in the coupling 
the form: 



and 



g 2 e 2 

'Afn 



6ag 2 (2j + l)e 2 
hE 2 ^ 2 



+ 0(a 2 ), (25) 



( ■ 



(27) 
\ 



+ jsi + ±aj(j + i) + w x (j, n) - Wo (j, n) A(j) 

A(j) + (j + 1)0 + J-a(j + l)(j + 2) + Witf + 1,0) + Wb(3 + 2, O) 



for the states |je) and |(j+l)g). The ground state |0) e ff = 
|0g) is an eigenstate of Ti e s with eigenenergy: 

E = ft(-A 6 /2 + Wi(0,n) + W (l,fi)). (28) 

Due to the doublet structure the blocks of the effective 
Hamiltonian are 2x2 matrices and the corresponding 
eigenvectors are for j > 0: 

|2j+l) eff = cos (|) | (j + 1) 5 ) + sin (|) lie), (29) 

|2j + 2) efi = -sin(|)|a + l).g>+cos(|)|je), 

where tanrfo- = 2 ^ J ^ and < rjj < n. Moreover, 

Ab _ n - 3a{ i + 1} + WtU, si) - w x (j + i, n) 



-w (j,ci)-w (j + 2,n) 



(30) 



In turn the eigenstates of the qubit-nonlinear oscillator 
system are obtained from the transformation 



\n) = exp(-iS)\n) eS . 
Finally, the cigcncncrgics are then 



(31) 



E2 j+ i/2 j+ 2 = Kj + \)n + \<*{j + 1) 2 

+h(w 1 (. h n) + w 1 (j + i,n))/2 

-hW (j,n)/2 + hW {j + 2,fl)/2 



+ -^| + 4|A(j)p. 



(32) 



These eigenergies are also eigenenergies of Htls-no by 
construction and are depicted in figure [T] (solid lines) for 
the case of an unbiased TLS, e = 0. At finite coupling 
the degeneracy is lifted and we observe avoided crossings 
(solid lines in figuref!]) . Due to the coupling the resonance 
condition acquires a shift compared to p^]) , the so-called 
Bloch-Siegert shifts, 



n = 



Ab ~ h a[j + + Wl{j ' Afc) ~ Wl{] + x ' Afc) 



+3 



W Q (j,A b )~W (j + 2,A b ) 
^9 2 A 2 



2hA\ 



ij + iy + 0(a 2 ig *). 



(33) 



The resonance corresponds to Sj = 0. We notice that the 
effect of the nonlinearity onto the Bloch-Siegert shift is 
very weak, namely at least of order 0(ag 2 ) and negligible 
for the values of nonlinearity and coupling we considered 
in the following. 

At resonance, equation (|33p . the minimal splitting of the 
former degenerate gap is: 



E 2 j+2 - E 2j+ i = ht/j + lg 
+0(a 2 ,g 3 



A 



a(j + 1) 



(34) 



We notice that at any point of our calculation we can 
set the nonlinearity to zero and reproduce the results 
obtained for the TLS-linear oscillator system^. 
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B. Dynamics of the qubit for the non-dissipative 
case 

The time evolution of the qubit-nonlinear- 
oscillator system without bath is given by 
pit) = exp(-|7Y T LS-No)p(0)exp(+^7i T LS-No) and 
therefore 

Pnm(t) = (n\p(t)\m) = exp(—iw nm t)p nm (0), (35) 

where u nm = i (E n ~ E m ). Consequently we obtain for 
the population difference in 



Pit) = po+J2 

Pnm (0) cos u nm t, 



(36) 



where we introduced po = Y^ n Pnn(0)- We observe from 
([56]) that the dynamics of the TLS is determined by an in- 
finite number of oscillation frequencies rather than show- 
ing a single Rabi oscillation. To set the initial conditions 
we assume that the qubit starts in the state \R) and that 
the occupation numbers of the NO are Boltzmann dis- 
tributed: 



p(0) = \R)(R\ 



1 



2no 



exp(-/?H NO ), 



where 



Zno = 5>x P [-/?(ftnj + 2 Q J'(i + 1))] 

3=0 



(37) 



(38) 



is the partition function of the oscillator and (3 = 
(fcsT) -1 is the inverse temperature. In the TLS-NO 
eigenbasis we get: 



Pnm (0) 



(n\p(0)\m) (39) 

1 00 3 
— exp[-/3(nfij + -aj(j + 1))] 

ms 1 7T ) Mjs) + sin (?) Hj' e 



2 

cos ( - ) /»> 



2 , 

sin ( ^ ) O'ei"' 



1, Low temperature approximation 



Equation (|36|) allows us to describe the non-dissipative 
dynamics in terms of the approximate eigenenergies and 
eigenstates (|3"Tj) and (|32l) . which involve in this way all 
nonlinear oscillator states. Therefore the Hilbert space 
under consideration is infinite. To calculate Pnm(0) and 
Pnn{0) we need to know the structure of a matrix el- 
ement such as (j,{g/e}\n) = (J, {g/e}| exp(-iS)|n) efI . 
The |n) e fi are themselves linear combinations of the un- 
coupled states \ j, {g/e}), see (|2"9")1 . Because we calculated 
exp(— iS) up to second order in the coupling Hamilto- 
nian 7ii n t, we find that the oscillator index j can at most 



change by four, see appendix [XJ For typical experiments 
on qubits the temperature is restricted to the regime of 
/3 _1 <C Ml, HAb- Due to the exponential function in (j3U)) 
high levels of the NO are only weakly populated and con- 
sequently we can truncate the infinite sum in equation 
(|39|l for the matrix elements of the density matrix at ini- 
tial time to j = 1. This means that the lowest 12 {|n}} 
states enter (|3T))) . 

After a close analysis we observe, by inserting (|39|) into 
(fTS"|) . that the coefficients p nrn (0) with n > 7 are of higher 
than second order in g. The same is valid for P50j£>60)P55 
and P66- Thus those terms do not occur in the calcu- 
lation of Pit). Of the remaining contributions we ob- 
serve that those with n — 5, 6 are either at least of order 
g exp[— /J(M2 + 3a)] or of order g 2 exp[— /3(7tfi + 3a)] or 
of order ag 2 . Thus we can also disregard contributions 
from p nm for n > 5 for the parameters chosen in the 
following, i.e., in the considered low temperature regime 
it is enough to restrict to the five lowest eigenstates of 
Wtls-no- Therefore the number of possible oscillation 
frequencies w nm is reduced to 10, where n, m = 0, 1, . . . , 4 
and n > m. 

In the following we show the dynamics of an unbiased 
TLS (e = 0), which results in vanishing of po, P3o(0), 
P4o(0), P2i(0) and £43(0). Therefore we obtain: 



Pit) = piocos(wiot) + P20 cos(w 20 i) 
+p 3 i cos(w 3 ii) +p 4 i COS(W4it) 
+P32 C0s(w 32 i) + P42 COs(w 4 2i)- 



(40) 



Exemplarily we consider in the following the resonant 
case for the corresponding linear oscillator, where f2 = 
A& = A . This corresponds to a slightly detuned 
nonlinear-oscillator system. The resulting transition fre- 
quencies using (j32j) are: 



W10 

W20 
W31 

W32 
OJ42 



n- 9 



3a 
2h 



3a 



9ag 9ag 2 

mi ^ 

9ag 



Ann 2 '- 

9ag 2 

4nn 2 '' 



n + ,g(l - 
n + gil + 
Q - gil + 
Q - gil - V2) 



2) + ^ 



9a 


9ag ' 


2h + 


4M1 l 


9a 


9ag ' 


2h ~ 


AMI V 


9a 


9ag ' 


2h + 


4M1 I 


9a 


9ag ' 


2H ~ 


AMI V 



2V2- 1 
2V2 + 1 
2V2 + 1 
2V2-1 



(41) 



9ag 2 
2htt 21 

2hn 2 ' 

9a 5 2 

2htt 21 

Qag 2 

2hn 2 ' 



Due to the nonlinearity the six different oscillation fre- 
quencies in equation (|41[) are shifted to higher frequencies 
compared to the linear oscillator case a = 0. In contrast 
to the linear case they are no longer located symmetri- 
cally around O = Aq. The reason for this lies in the 
non-equidistant energy levels of the nonlinear oscillator 
alone and in the interplay of coupling and nonlinearity. 
The population difference P(t) and its Fourier transform 
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are shown in figure [2] As in the linear case, the dom- 
inating frequencies are uj w and to 2 o- These correspond 
to transitions between the first and the second state of 
the qubit-NO-system and the ground state. In the linear 
oscillator case the weight of their peaks is almost equal, 
whereas with weak nonlinearities the peak correspond- 
ing to wio is more pronounced. This is due to the fact 
that the frequency corresponding to the more pronounced 
peak fits more accurately the resonance condition, which 
includes the Bloch-Siegert shift in (J33J). The weight of 
the peaks can additionally be influenced by allowing a 
finite bias of the qubit, e / 0. The zero bias case was 
chosen here for simplicity. 

From these graphs and equations (f34|) and (|4Tj) we can 
read off first that the vacuum Rabi splitting is decreased 
for finite nonlinearity and second that the overall fre- 
quency shifts compared to the linear case are larger the 
higher the oscillator levels are involved if the coupling 
g is not too large to overcome the effects caused by the 
nonlinearity. 



IV. INFLUENCE OF THE ENVIRONMENT 

The knowledge about decoherence and dissipation pro- 
cesses entering in the qubit dynamics is essential for 
quantum computation. Therefore we consider now the 
qubit-nonlinear-oscillator system to be coupled to an en- 
vironment and treat the full Hamiltonian TL. 

A. Master equation for the qubit-NO system 



As shown in section III BL equation (fT4")) , we need for 
the calculation of P(t) the density matrix p(t) of the 
qubit-nonlinear oscillator system. To take into account 
the effect of the bath we start from the Liouville equation 
for the full density matrix W(t) of Tt, 



lh -dT 



[H 



NO-BJ 



(42) 



where the index / denotes the interaction picture. 
Following 5 —^ we arrive at a Born-Markov master equa- 
tion for p(t) being in the Schrodinger picture and ex- 
pressed in the basis of the eigenstates of 7Yq-no : 

(t) + C nm ,kiPki(t)- (43) 

k,l 

The first term includes the free dynamics, whereas the 
second accounts for the dissipative one. The Bloch- 
Rcdfield tensors are defined by: 

= [G(to nk )N nk - G(u)i m )N m i] (44) 
v 

lyik'Vk'raj 
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FIG. 2: (Color online) Top: Dynamics of the population dif- 
ference P(t) for the unbiased, e = 0, qubit-nonlinear oscilla- 
tor system at linear resonance (fi = Ao) (blue (dark gray) 
line). We choose a nonlinearity a = 0.02M1, a TLS-NO cou- 
pling g = 0.18f2, and inverse temperature f3 = lQ(Kt)~ ■ For 
comparison we plotted the corresponding linear case (orange 
(light gray) line). Bottom: Fourier transform F(uj) of P(t) 
for the unbiased system. The dominating frequencies are wio 
and ui20- To visualize the delta- functions, finite widths have 
artificially been introduced. 



with N nm = | [coth(ft/3w nm /2) - 1] and y nm = (n\(B + 
B')\m). In the following we assume to have an 
Ohmic bath described by the spectral density G(ui) = 

GohmM = KU>. 

For the derivation of the master equation besides the 
Born-Markov approximation more assumptions have 
been made. We only mention them briefly: first, we 
assume that the system and bath are initially uncor- 
rected (at t = 0), i.e., W(0) = pi(0)p B (0), where 
Pb(0) = Z^ 1 exp(— /37Yb) and Zq is the partition func- 
tion of the bath. Because the bath consists of infinite 
degrees of freedom we assume the effects of the interac- 
tion with the TLS-NO system on the bath to dissipate 
away quickly, such that the bath remains in thermal equi- 
librium for all times t: W\(t) = pi(i)ps(0). Additionally 
an initial slip term is neglected, which occurs due to the 
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sudden coupling of the system to the bathiii. Finally we 
disregarded the Lamb-shift of the oscillation frequencies 

B. Matrix elements 

The Redfield tensors, equation (|44]) , depend on the 
matrix elements y nm of the NO position operator in the 
TLS-NO eigenbasis. Using equation (|3"Tj) we rewrite y nm 
in the form: 

Vnm = (n\y\m) = e s(n\ex.p(iS)yexp(-iS)\m) eB 

= c tf(n\y\m) cS . (45) 

The effective states are given in (|29|) as linear combina- 
tions of states of the {|jg); |je)} basis. In the following 
we show the different building blocks for y nm . We can 
distinguish between different situations. First there are 
matrix elements where neither the qubit nor the oscilla- 
tor state is changed, namely: 

(jg\y\jg) = -2(L L Oo(g) + L NO0 (j,at,g)), (46) 
(je|y|je) = +2(L LO o{g) + L NO0 (j,a,g)), 

where L LO o(g) = gs/A b tt and L NO a{j,a,g) = 
—6age(2j + l)/HAb^ 2 . These matrix elements contain 
contributions independent of the oscillator occupation 
number j for zeroth order in the nonlincarity a and ac- 
quire a level dependence in first order. 
A transition within the qubit is described by 

(jg|y|ie) = L LO0+ (g)+L NOQ+ (a,g){2j + 1). (47) 

Here we introduced abbreviations, given in appendix 151 
to show the actual order of the matrix elements involved. 
The notation is as follows: indices LO and NO refer to 
the linear or nonlinear oscillator, respectively. An ad- 
ditional index number, Aj, indicates that the nonlinear 
oscillator state is changed by Aj quanta. We have ele- 
ments where zero, one, two or three quanta are emitted 
or absorbed by the oscillator. Moreover we introduce 
indices +/— or g/e which correspond to the TLS transi- 
tion g — > e or to e — ► g, respectively, or to the qubit not 
changing from g or e configuration. 
For the case Aj = 1: 

(ig|y|(i + i)g) = V7+T[i + (j + i)L NO ( a )+ 

Lloi{9 2 ) + L NO i g {j,a,g 2 )] , 
(My\(j + l)e) = y/j + 1 [1 + {j + l)L NO (a)- 

L L oi{g 2 ) + Ljvoic(j,a,g 2 )] ,(48) 

(jg|y|(i + i)e) = ^JTT[L L01+ (g 2 )+ 

L NO i+(a,g 2 )(j + 1)] , 

(je|y|(j + l)g) = y/j + T[L LO i-(g 2 )+ 

L NO i-(a,g 2 )(j + l)] , (49) 



describe processes where an oscillator quantum is ab- 
sorbed. All the matrix elements in (|47[) . (|4"5|) and in 
(|49|) contain both zeroth-order as well as first-order con- 
tributions in the nonlinearity. Additionally, due to the 
fact that the states of the NO are linear combinations of 
the linear oscillator states, see equation ([7]), additional 
transitions involving a change of the oscillator state by 
more than one quantum are allowed. They correspond 
to Aj = 2, Aj = 3 and read as 



O'g 


;|y|CH 


-2)g) 


= vu- 


M)(i- 


- 2)L N02 (a,g), 




>\y\ti - 


f 2)e) 


= Vti- 


i- 1)0' - 


- 2)L N0 2+(a,g) 



(je|y|(j + 2)g) = y/(j + + 2)L N02 -(a,g), 
(je\y\(j + 2)e) - - y/(J + l)(j + 2)i w02 (a, g), 



0'g|y|(j+3)g) = V0' + l)0" + 2)0' + 3) [L N03 (a,g 2 ) 
-L NO {a)/2] , 

O'g|y|0'+3)e) = y/(j + l)(j + 2)(j + 3)L JV03+ (a, g 2 ), 

(je|y|(j + 3)g) = y/(j + l)(j + 2)(j + 3)L N03 -(a, g 2 ), 

0'e|y|(i + 3)e> = V(j + l)(j + 2)(j + 3) [-L N03 (a, g 2 ) 
-L NO (a)/2] . 

Notice that all terms in (T5D|) vanish when a — 0. The 
terms in (|48p and (|50|) involving no change in the qubit 
and a change in the oscillator by Aj = 1 and Aj = 3 
quanta contain g- independent nonlinear contributions. 
The interplay of nonlinearity and coupling in lowest order 
can be observed in (jg|j/|je), and in the terms involving 
an oscillator level change by 2. Additionally at the degen- 
eracy point, e = 0, L LO o(g), L NO o(j,a,g), L LO i±(g 2 ), 
L N0 2(a,g), L NO i±(a,g 2 ), L N0 3±(a,g 2 ), and parts of 
Ljsroi{ g /e}{j, a, g 2 ) vanish. We are now able to calculate 
the matrix elements y nm - They are given in appendix IB1 



C. Dissipative dynamics 

To calculate Pit) we have to solve the system of cou- 
pled differential equations (|43]) . When several TLS-NO 
levels are involved an exact solution can only be found 
numerically. Hence, in the remaining of this section we 
discuss three different approximation schemes, two based 
on the full secular approximation (FSA) applied to ((43|) 
and one based on a partial secular approximation (PSA). 
We then compare the so obtained analytical predictions 
with the exact numerical solution of (|4"3"1) . 



1. Full secular approximation (FSA) 
We define: 

Pnm(t) = exp(-iuj nm t)a nm (t), (51) 
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which, inserted in (|43|) . enables us to obtain a set of dif- 
ferential equations for a nm (t): 

o"nm(*) = n^2£nm,kl exp[i(w nm - u k i)t]a M (t). (52) 

kl 

The FSA consists of neglecting fast rotating terms in 
equation ([52jl such that only terms survive where u nm — 
LJki vanishes. This allows an effective decoupling of diag- 
onal and off-diagonal elements such that 

kk<Jkk(t), (53a) 

k 

Onm{t) = Tr£nrn,nmO-nm(t) for n ^ 1X1. (53b) 

The off-diagonal elements are determined by 

(?nm(t) = cr° m exp(7r£„„ l! „„ l t), (54) 

which results with ([BT]) in 

Pnm(t) = p° nm exp(ir£ nm , nm t)exp(-iuj nm t). (55) 

The separation of the oscillatory motion of the dynamics 
from the relaxation one allows us to divide (fT4|) into two 
parts 

P(t) = PrelaxW+^d cphas (t), (56) 

where P re iax(0 = ^2 n Pnn{~k) is the relaxation contribu- 
tion and Pdcphas(i) = 2 n > m Pnm(t) is tne dephasing 
part. Inserting (|55p in the last expression and using (|15p . 
we obtain: 

-fclcphas(^) — ^ Pnm 

(0) exp(-r t)cos(uj nm t), (57) 

where the dephasing rates are determined by T nm = 
—itC nm nm . The actual form of the dephasing coefficients 
£nm,nm can be found in appendix [Cl and the initial con- 
ditions /9„ m = <r° m = p nm (0) are defined in (f3"9")) . The 
diagonal elements are more difficult to obtain, since the 
coupled system of differential equations in (|53a[) has to 
be solved. To proceed we restrict ourselves in this sec- 
tion again to the physical relevant low temperature case, 
such that the highest qubit-nonlinear oscillator state in- 
volved is the eigenstate |4). Calculating the rate coeffi- 
cients accompanied with the five lowest eigenstates, we 



observe that there are only eight independent ones due 
to the structure of the rate coefficients. These are £qo ii> 
£00,22, £11,22, £11,33, £11,44, £22,33, £22,44, and £33,44- 
In general they are determined by 

£jj,fcjfc = %G(uj : jk)N ] ky 2 j k with j < k, (58) 

where j and k adopt the above values. Furthermore, 
£00,33, £00,44, £33,00 and £44,00 are disregarded, because 
they are at least of order 0(g i ). The remaining rate 
coefficients are combinations of the above. We find that 

£/cfe,jj = £jj,fefe + 2G(cj J - /c )?/| fe (59) 
= (N jk + l)2G(u> jk )y%, 

and 



£00,00 = 


— £11,00 


— £22,00 






£11,11 = 


— £00,11 


- £22,11 


— £33,11 


— £44,11 


£22,22 = 


— £00,22 


— £11,22 


— £33,22 


— £44,22 


£33,33 = 


— £11,33 


— £22,33 


— £44,33 




£44,44 = 


—£11,44 


— £22,44 


- £33,44 





Low temperature approximation (LTA) 

Despite the above relations equation (|53a[) is too com- 
plicated to be solved analytically. Therefore an addi- 
tional approximation is applied: we consider the fac- 
tor N nm + 1 = | [coth(hf3uj nm /2) + 1] with n < m in 
equation (|59"|) and use that lim a ,^_ 00 coth(a;/2) = —1 is 
reached exponentially fast. 

The terms containing this factor are neglected in the fol- 
lowing. As we consider only the lowest five levels, this 
amounts to require max{w„ m } = |wi4| 3> k^T. Using 
equation (|34[) we observe that W12 oc g and W34 cx g. For 
this reason and due to the structure of y nm given in equa- 
tion (|Bip the rates £11,22 and £33,44 are at least of order 
0(g 3 ) and can be neglected. With equation the rate 
matrix £ rc iax associated to (|53a|) becomes 



(0 


£00,11 


£00,22 





\ 





- £00,11 





£11,33 


£11,44 








- £00,22 


£22,33 


£22,44 











—£11,33 — £22,33 





\o 











—£11,44 — £22,44 / 



r 



The eigenvalues and eigenvectors of this matrix and the in appendix [D] In contrast to the simple analytic expres- 
associated time evolution of the elements a nn (t) are given 
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sion for the dephasing part the relaxation rate is not easy 
to extract as P re i a x(i) = J2 n Pnn(t) consists of a sum of 
several exponential functions, cf. (fl5|) and appendix |D] 
However an analytical formula for P{t) can be provided 
using (|56j) . 



Smallest eigenvalue approximation (SEA) 



0.15 



< 

O , 
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FIG. 3: (Color online) The relaxation rate IV given in equa- 
tion (|62l) drawn against the oscillator fequency £1 (contin- 
uous blue (dark gray) line). We used e = 0.5Ao, corre- 
sponding to a frequency splitting At = I.II8A0, coupling 
g — O.I8A0 and the nonlinearity a = 0.02ftAo. The damping 
constant is k, — 0.0154 and /3 = 10(SAo) _ . At resonance 
2 A 2 

{Q = A b - 3a/ h + 2^0.) r r is maximal. For comparison also 

the second lowest eigenvalue is plotted (orange (light gray) 
dashed line). The inset shows the two eigenvalues close to 
resonance. 



In order to get a better insight into the effect of the 
relaxation mechanism, we consider the long-time dynam- 
ics of the system. This means that the we direct our 
attention to the smallest eigenvalue of the relaxation co- 
efficients, which dominates at long time, rather than to 
tackle the many relaxation contributions involved in the 
populations a nn it) . We do not make the low temperature 
approximation discussed above. We restrict for simplic- 
ity to the three lowest qubit-NO eigenstates |0), |1), |2) 
in (I53aj) and obtain using ((60|) : 



(61) 



£11,00 — £22,00 


£00,11 


£00,22 


£11,00 


—£00,11 — £22,11 


£11,22 


£22,00 


£22,11 


—£00,22 — £11,22 



We do not neglect £11,22 and £22,11, even if they are 
at least of order 0(g 3 ), because these contributions lift 
the degeneracy of the two lowest eigenvalues at resonance 



(see figure [3]). The smallest eigenvalue is 

r r = 

£ 



(62) 



^ n .in in 



— 4(£oO,ll£oO,22 + £ll,Oo£oO,22 + £oO,ll£ll,22 
£ll,Oo£ll,22 + £(X),ll£22,00 + £ll,22£22,00 

n l/2 

+£22,ll£oO,22 + £ll,Oo£22,ll + £22,0q£22,11 ) 



Additional detuning allows for a further simplification: 

2 A 2 

r r sa 7r£oo 22 for fl + 3a/h~ < Ab and IV as 7r£ o 11 

a 2 A 2 

for Q + 3a/h— A3 ° > A;,. In figurc[3]thc relaxation rate 

r r in ((nil) is plotted as a function of the linear oscilla- 
tor frequency Q. It is maximal at resonance, whereas it 
decays for Q being detuned from resonance. Addition- 
ally we plotted the second smallest eigenvalue of (fBTj) for 
comparison (dashed orange (light gray) line in figure [3]). 
In the long-time limit it then holds: 



frelax(*) = (j?0 ~ Poc)e FA + Poo, 



(63) 



where, like in section UlI B[ po = ^2 n Pnni0). To obtain 
Poo we have in principle to find the steady-state solution 
of (|53ap . Here, we just assume for t — > 00 a Boltzmann 
distribution for the TLS-NO system, so that p n „(oo) = 

^TLS-NO ex P( _ ' 9£, «) witn ^TLS-NO = J2n ex P(-/^n)- 

Thus, 



(jsH 2 - (je\n)' 



(64) 



Poo = \ cos e 

» 3 

+2 sin 6 ( jg I n) (j e| n) } p nn (00) . 

The formula for the long-time dynamics is then obtained, 

Pit) = (j5 -Poo)exp(-r r t)+p oo + (65) 
^ jw(0)exp(-r nm i) cos{u nm t). 



To get further insight on the dominant frequencies we 
evaluate the Fourier transform of (|fJ5|) according to 



/>oo 

F{oj) = 2 dt cos ujtP 'it), 
Jo 



(66) 



yielding 

F(oj) = 
2(po ~Poo) 



(67) 



^ 2 + r2 
1 



(o)r 
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2. Partial secular approximation (PSA) 

The PSA is an improvement to the FSA, where one ac- 
counts for corrections to the equations for the coherences 
due to dominant rotating terms ui nm — u^i in (|52p . The 
equation for the populations is still given by (|53aP ■ At low 
temperatures the dominant correction to the FSA comes 
from transitions involving the quasi-degenerate states |1) 
and |2). To solve the off-diagonal part we have to de- 
termine (Toi, (702; 013, 023, <T 14; an d &24- With ([52]) the 
system of equations is 

Pnm(^) — ( i^nrn £>nrn.nm) Pnmij^) ~t~ £*nm,jkPjk{^) , 

(68) 

I 



We can calculate analytically the relaxation and dephas- 
ing part of P(t). While the FSA allows a simple form for 
the dephasing rates, T nm = — Tr£ nm nm , the PSA one is 
much more involved. As in case of the SEA the smallest 
eigenvalue dominates the dephasing behavior. The cor- 
responding Bloch-Redfield tensors are found in appendix 

o 



V. NUMERICAL VERSUS ANALYTICAL 
PREDICTIONS 

In the following we compare the results for the dynam- 
ical quantity P(t) and its Fourier transform, obtained by 
a numerical solution of (|43[) , with the predictions of the 
approximations from section ITVl 



A. Low temperature 



Pjk(t) = nCjk, nm Pnm(t) + {—iujk + n£jk,jk)Pjk{t) (69) 

with {(nm),(jk)} = {(01); (02)},{(13); (23)}, or 
{(14); (24)}. The solution is 

Pnm — C nm,jk V nm,jk eX P(^nm,jk^) (?0) 
+ c nm,jk v nm,jk ex P(^L2,jfc*)> 
Pjk = c nm,jk ex P(^im,jfc*) + C nm,jk ex P(^Ln,jfc*)) 

where the oscillation frequencies and the decay of the 
off-diagonal elements are given by^ 



(71) 



(72) 



(73) 



(74) 



LTA, and PSA) to the numerical solution in figure [4] 
We recognize that the dynamics and the corresponding 
Fourier spectrum are well reproduced within the simple 
SEA approach as well as in the two LTA and PSA treat- 
ments and determined by the superposition of two os- 
cillations. The best approximation is the PSA. In the 
following we use the SEA approach due to its simpler 
analytic form. 

To determine the effects of the nonlinearity onto the qubit 
dynamics we compare P(t) and F(u>) with the corre- 
sponding linear case in figure[5] We choose Q = A&. Both 
in the nonlinear and in the corresponding linear case two 
oscillation frequencies are dominant. Due to the Bloch- 
Siegert shift, see (|3"3"|) . in both cases = is not the ex- 
act resonance condition. However, in the nonlinear case 
the nonlinearity partly compensates the Bloch-Siegert 
shift, which also influences the relative peak heights, as 
we argued in section IIII B[ 



with 

Pnrnjk \J [^{^nm.nm £jk,jk) ^(^nm ^j/c)] ~t~ ^ 1 ^'^^~ , 7unjk^~-jk,nm- 

The amplitudes of the oscillations are given through the coefficients 

(+/ — ) .^^^jk^nmPyun Pjk [^(.^nm^nm £jk.jk) ^(^nm ^jk) "F Pnmjk] 

^nrnjk n ry 

^-^nmjk 

and 

V nm,jk = 9 r M-Cnm.nm _ £jk,jk) ~ i{^nrn — OJjk) ± Rnm,jk] ■ 



We start by focusing on low temperatures j3 — 10/ (Ml) 
and compare the results for all three approaches (SEA, 
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FIG. 4: (Color online) Comparison of the behaviour of P(t) 
and its Fourier transform F(ui) as obtained from the numer- 
ically exact solution (red (dark gray) line) and the three ap- 
proximation schemes (orange (light gray) line) discussed in 
the text. Top: Smallest eigenvalue approximation (SEA), 
Middle: Low temperature approximation (LTA) and Bottom: 
Partial secular approximation (PSA) . The chosen parameters 
are: a = 0.02M7, g = 0.18O, e = Ofi, k = 0.0154 and 
j3 = 10(M1) _1 . The dynamics is well reproduced within all 
approximations, however the agreement of the PSA with the 
exact numerics is the best. In the corresponding Fourier spec- 
trum almost no deviations occur for all three approaches. 
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FIG. 5: (Color online) Top: P(t) for the linear (orange (light 
gray) line) and nonlinear (blue (dark gray) line) case using 
the parameters: a = or a = 0.02/ifi, respectively, and g = 
0.180, k = 0.0154, e = 0, A b = tl, (3 = 10(M)~\ Bottom: 
Corresponding Fourier transform F(uj). 



B. Higher temperatures 

To investigate the influence of temperature we show in 
figure [S] P(t) and the corresponding F(oj) for the same 
parameters as in figure [5l but at inverse temperature: 
f3 = 3/(hfl). By increasing the temperature higher oscil- 
lator levels are populated and influence the dynamics of 
the qubit. We calculated the corresponding equations for 
the long time dynamics within the SEA. The relaxation 
matrix jC re iax for the rate T r was calculated inplementing 
higher levels, until |8) c ff- 

We plot for comparison also the linear oscillator case. We 
observe again the overall shift of the resonance frequen- 
cies to higher values and that a new shoulder arises in 
the Fourier spectrum. It corresponds to the transition 
frequency uj 32 (see also figure [5] bottom) . We checked 
numerically that the structure of the Fourier spectrum 
can be fully respresented by summation of the six contri- 
butions in P{t) with the frequencies wio, W20, ^32, W42, 
W13, and W14. These six contributions arise due to the 
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finite populations of the involved levels. Therefore the 
appearance of additional shoulders in the dynamics is a 
pure temperature effect, which is also seen in the cor- 
responding linear case. The frequency shift induced by 
the nonlinearity is much larger for the higher levels. The 
effect of temperature is also reflected in the height of the 
dominating peaks, which is decreased for higher temper- 
atures. The temperature can not influence which peak is 
dominant. This means by comparing figure [5] with figure 
[5] that in both figures in the nonlinear case the peak cor- 
responding to u>iq dominates over the one corresponding 

^20- 

The use of a nonlinear oscillator instead of a linear one 
has advantages which rely in the fact that the energy 
spectrum of the nonlinear oscillator is not equidistant. 
Supposing that the TLS frequency A& can be tuned, it 
is in case of the nonlinear oscillator possible to have the 
TLS in resonance with exactly one and only one nonlin- 
ear oscillator state transition. All other transitions are 
then off resonance/detuned. For the linear oscillator in 
resonance with the TLS the number of possible transi- 
tions is in principle infinite. Therefore we determine in 
the following the dynamics of the qubit by putting the 
qubit in resonance with the nonlinear oscillator transition 
1 3) — > 1 2), see figure [71 We read off from figure [7J that 
the detuning compared to figure [5] results in the enhance- 
ment of the W2o-peak, whereas the other dominating peak 
is shrinked. This is due to the different resonance con- 
ditions leading to opposite weights of the peaks for the 
nonlinear case in figure [7] compared to figure [5J However 
a peak corresponding to higher transitions is not seen. 
The reason for this is the small population of the higher 
oscillator levels involved. 
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FIG. 6: (Color online) P(t) and its Fourier transform F(ui) for 
the parameters: a = 0.02M7, g = 0.18S1, « = 0.0154, e = 0, 
A;, = Q as in figure but /3 = 3(h£l)~ . For comparison we 
plotted the linear case in orange (light gray). 



VI. CONCLUSIONS 



To conclude, we determined the dynamics of a TLS 
which is coupled via a nonlinear oscillator to an envi- 
ronment described by an Ohmic spectral density. We 
restricted ourselves to the regime of weak nonlinear- 
ity weak damping and moderate coupling of oscillator 
and TLS. To diagonalize the qubit-nonlincar-oscillator 
Hamiltonian we used Van Vleck perturbation theory, 
hence avoiding the use of the rotating wave approxima- 
tion (RWA). Within the RWA and for vanishing nonlin- 
earity our model would reduce to the Jaynes-Cummings 
Hamiltonian. In section IIII Bl an analytical expression 
for the non-dissipative dynamics is given, which accounts 
for the infinite Hilbert space of the composed system. 
The influence of the nonlinearity onto the qubit dynam- 
ics is determined and compared to the linear case. 

At low temperatures ksT < M7, ftA/, this infinite 
Hilbert space can be truncated such that the transition 
processes between the ground state and the two first ex- 
cited energy levels of the qubit-nonlinear-oscillator sys- 
tem dominate the dynamics. As in the linear case this 
yields a pronounced vacuum Rabi splitting. 



To investigate the influence of the bath we solved 
the Bloch-Rcdfield Markovian master equation for the 
density-matrix of the qubit-nonlinear-oscillator system 
numerically. For an analytical treatment we considered 
three kinds of approximations: first a full secular ap- 
proximation including a low temperature approximation, 
where all fast oscillating terms are neglected. Second an 
ansatz for the long-time dynamics allows a general ex- 
pression for the relaxation and dephasing rates of the 
qubit. The third approximation was a partial secular ap- 
proximation reproducing almost perfectly the exact nu- 
merical solution. A comparison of these three analytical 
approaches showed good agreement with the numerical 
solution. Finally, we investigated the effect of the non- 
equidistant energy spectrum of the nonlinear oscillator 
on the TLS dynamics. To do so, we allowed higher tem- 
peratures to populate higher levels and moreover we con- 
centrated on the actual transition of the nonlinear oscil- 
lator from 1 3) — > 1 2). We observed the rise of additional 
shoulders in the Fourier spectrum and showed that the 
shift in the transition frequencies is much larger if higher 
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200 



oscillator levels are involved. 
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FIG. 7: (Color online) P(t) and its Fourier transform F(u) 
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APPENDIX A: VAN VLECK PERTURBATION THEORY 



In our case the perturbation Hj nt is proportional to B + . Therefore we consider first the action of this operator 
on arbitrary nonlinear oscillator states \l), \ m): 



{l\B + B^\m) = (l\ y/m\m - 1) + a ( 2 m) Vm + 2\m + 1) + a%Wm - 2\m - 3> + a [m 4 Wm - 1\m - 5) (Al) 
+4 m Vm + 4|m + 3) + Vm + l\m + 1) + 4™Vm + 3|m + 3) + a^Vm- l\m - 1) 
+a!_™Vm- 3 1 m - 3) + a^y/m + ^m + 5) 1 + C(a 2 ), 
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where \l)o denotes an eigenstates of the linear oscillator. Now we have different cases: 



I = m — 1 

I = to- 3 

I = to — 5 

I = m+1 

I = to + 3 

I = TO + 5 



(m- 1|(B + Bt)| 
(to — rn 
{m-5\(B + B^)\m : 
(m + l\(B + B r )\m] 
(m + 3\(B + B r )\m 



= V^+aL7V^^ + 4 ) V^Tl + 0(a 2 ), (A2) 
= a ( ™ ) V^2 + a^V^Z + 4 m ~ 3) \/^ + ai m - 3) V^TT + 0(a 2 ), 
= + 4 m - 5) + O (a- 



0, 



- V^TT + 4™ +1) \/to" + 4 m) V^T2 + G(a 2 ) 



a (m + 3) V^+T +a ( Z +3) V^+a ( r ] + + O (a 2 ) , 

(m + 5|(B + Bt)|m) = 4™ +5) V^TT + 4 m) V^T5 + G(a 2 ) = 0. 



Due to the manifold structure we only have to consider for Van Vleck perturbation theory the matrix elements 
involving Z = to±1,Z = to±3. Therefore we introduce the notations: 



VJ+T 



2M] 



a(j + 1) 



(A3) 



ns(j,a) = a« + a« V7^3 + 4 j " 3) V7 + a£' _3) VJ + T = - 1)0' - 2). (A4) 

The non-vanishing matrix elements for the transformation matrix are in first order: 



,-o(i) 

^(j-l)e,je 



(e,j- l|Hin t |e,j) = g^^i(i-l) = geyj 

£ e (j-l) - £ ej + ^-a • 2j A b n 



1 " 27^' 



+ 0(a 2 ), 



(A5) 



(i) 



(g,i|Wi„ t |g,j + i) 



ig,0'+i)g R . _ R . . , 
-^gj ^gO+i) 



A a . 2 (j + 1) 



+ G(a 2 ), 



(l) = (g, j|Hi nt |e, j + 1) 

^gj ^cO + l) 



5st«i(j) 



= gAoVJ+T 

A 6 + 0+^a-2(j + l) A b (A b + fi) 



3a(jJ-l](Afe + 30) 

2M}(A b + 



0(a 2 ), 



(e,j|Wi nt |e,i + 3) _ 5^n 3 (i + 3,o) + ^ 2) 



je,(j+3)e TP . _ j? . 



317 



,-o(i) 

io ig,(i+3)g 



(g,j|Wint|g,j + 3) ^n 3 (i + 3,a) 



E gj - E, 



g(j+3) 



3ft 



+ G(a 2 ), 



^g,0+3)e 



(g,j|Wint|e,j + 3) _ 3^3(i + 3,a) 



^gj _ E e(j+3) 



A b + 30 



+ 0(a 2 ), 



,-o(i) 

i,3 je,(j+3)g 



(e, j|Hi„t|g, j + 3) _ Ssi^C? + 3, a) 



^ej - #g(j+3) 



-A 6 + 3fl 



+ 0(a 2 ). 



Due to the fact that 713 (j, a) is a purely nonlinear contribution, we can reduce the possible contributions for the second 
order of the transformation matrix. Restricting to lowest order in the nonlinearity the non-vanishing contributions 
are either combinations of involving twice n\{j) and expanding this afterwards to first order in the nonlinearity or 
combinations of both n\{j) and n^{j, a), while ni(j) is reduced in this case to the zeroth order in the nonlinearity 
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because nz(j, a) is already of first order in the nonlinearity. For the second order we obtain: 



lS jJ,(j+2)g 



E g(j+2) - E cj 



(e,j\Hi nt \e,j + 1) (c, j + l\H I nt \g,j + 2) (e,j\H i nt \g,j + l)(g,j + l|Wi nt |g, j + 2)" 



E e(j+1) - Eej 

(e,j\Hi n t\e,j-l)(e,j-l\H Int \g,j + 2} 

2 ( E g(j+2) - E ej) 
(c,j|^Int|g,j-l)(g,j-l|Hl nt |g,j + 2) 

2(E g(j+2 ) - E ej ) 
(g, J I^Int |g, j + 3) (g, j + 3|Hi nt |g, j + 2) 

2(£ s o+2) - E ej ) 
(e,j\Hi nt \e,j + 3)(e,j + 3\Hi nt \g,j + 2) 



1 



^eO-i) - E ej + E cU _ 1} - E g{j+2 ) 

1 1 

+ 



^g(j'-l) _ E *3 E g(j-l) ~ E g(j+2) 



+ 



2(25, 



(j+2) - ^cj 



E ej) 



E sU+3) ~ E ej E g(j+3) - E g (j+2) 



+ 



E c Q+3) - E ej E e(j+3) - E gU+2) 



2fiV £ Ao V /(j + l)(j + 2) 

A 2 b h 2 n(2n - A b ) 



3a(2j + 3)(A fc -30) 

nn(2n - A h ) 



g 2 (2j + 3)y/(j + l)(j + 2)aA e(A b - 50) 2 
12/10 2 (A* - 40Af + ft2 A 2 + 6Q3 Ab ) 



(2) 



Jg,0'+2)g 



(g, Jl^lntlg, j + l)(g, j + l|Hlnt|g, j + 2) 
2{ E g (j+2) - E gj) 

g,j|Wi nt |e, j + l)(e,j + l|Wint|g, j + 2) 



+ ■ 



£g(j+l) - E gj E g(j+1) - E g(j+2) 



( E g(j+2) - E g j)(E e ( j+1 j - Egj) 

g,j\n Int \ g , j - i)( g , j - i|Wi„ t |g, j + 2) 



2( E g (j+2) - Egj) 
g, j|Wint|g, 3 + 3)<g, j + 3|Wi„ t |g, j + 2) 



1 



2 (^g(j+2) - £gj) 

g, j|Wint|e, J + 3) (c, j + 3\Hm\g, 3 + 2) 



E gU-l) ~ E S3 + E g(j-1) - E g(j+2). 

1 1 

+ 



E sU+3) - E gj E g(j+3) - E g(j+2) . 



2( E g(j+2) - E gj) 

g,3\Hj nt \e, (j - l))(c,j - l\H lnt \g, (j + 2)) 



E c(j+3) - E gJ E c(j+3) - E g(j+2) 



{ E g(j+2) - E gj)( E c(j-l) - E g(j+2)) 



%V(J + !)(.? + 2) 



2A 2 o 



3ae 



A 2 , 



1 - 



+ 



_ 3a((2j + 3)A b + 0(3j + 4)) 

2h 2 n 2 1 h(A b + n) h(A b + n)n 

ffV(.7 + l)(j + 2)q / 2e 2 {{2j + 3)A 2 + 3(j - l)OA b - 3(j + 6)0 2 ) A 2 



8?iA 2 2 



+ 



(At -30) (A 2 + 40A b + 30 2 ) 



0(a 2 ), 
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*°je,(j+2)e 



(c,j\Hint\e,j + 1) (c, j + l\Hlnt\c,j + 2) 
2 (^eO+2) - E e j) 

(e, j|Hlnt|g, j + l)(g, j + l|Hint|e, j + 2) 

(E c (j+2) - E ej){ E g(j+l) - E e(j+2)) 

(c,j\Hint\g,j - l)(g, j - l|Hint|e, j + 2) 

2(-E , e (j+2) - E e j) 

(e, j|Hint|e, j - l)(e, j - l|7fr n t|c, j + 2) 

2 (E e (j +2 ) - E e j) 

(c, j^WK j + 3) (e, j + 3|Hi nt |e, j + 2) 

2(£'e(j+2) - E e j) 

<ej|Wi nt |g, j + 3)(g, j + 3|Wi nt |e,j + 2) 



+ 



E e (j+i) - E ej E e(j+1) - E c{j+2) 



Eg(j-l) - E cj E e (j-1) - E e {j+2) 



+ 



£ e(j-l) - E e j ^e(j-l) - E g(j+2) 



E c (j+3) - E gj E c(j+3) - E c (j+2) 



(E e (j+2) - E ej )(E g (j +S ) - E ej ) 



a 2 , 



ft(0 + A b 



1 _ 3a(A 6 (2j + 3) + n(3j + 5)) > j + 3ae 2 ' 



+ 



h{A b + Q)Q 

ff 2 aV(j + l)(i + 2) ( (3 + 2j)A 2 + 3(4 + j)A b Sl - 3(j - 3)r> 2 



8^A 2 2 



Ag + A 2 - m 2 A b - 90 3 



? -c(2) 



(gj j|Hlnt|g, j ~ l)(g, j - l|^Int|e, j) 

2 (-Eej - E g j) 

. (g>i|Wint|e,j + l)(e, j + l|Wint|e,j) 



1 



1 



E, 



2(E ej - Egj) 



s(j'-i) 
1 

E e (j+1) - Egj 



1 



E e 



O+i) 



E, 



+ 



(g,j|Wint|e, j - l)(e, j - l|Wint|e, j) (g, j|Wint|g, j + l)(g, j + l|Wint|e, j) 



(Eej ~ E gj)( E e(j-l) - E ej ) 



ftA 2 fi(A fc + n) 



h(2j + 1) 3a(2j 2 + 2j + l)(2A b + 30) 



(Ecj - Egj)(E g(j - +1 ) - E g j) 

+ 0(a 2 ), 



2n(A b + n) 



iS {2) 

ig,0'+2)e 



(gj Jl^lntlg, j + l)(g, j + l|^Int|e, j + 2) 

2(E ( i+2 ) - E gi ) 
(g,i|Wi nt |e,i + l)(e,i + l|Wi nt |e,j + 2) 



1 



+ 



- E g(j+1) ~ E B3 E g(j+1) ~ ^c(j+2) 
1 1 



2(E e (j +2 ) - E gJ ) 

<g,j|^Int|g,j- l)(g,j~ l|Hlnt|e, j + 2) 

2(E e( j +2 ) - Egj) 
(gj jj^intje, j + 3)(e, j + 3|7ft nt |e, j + 2) 
2(E e ( j+2 ) - E gi ) 



^(j+i) - Egj 



E e(j+1) - E c(j+2) 



1 



+ ■ 



S g(j-1) ~ E S3 E SU-1) ~ E e(j+2) 

1 1 

+ 



E e (j+3) - E gi 



E c (j+3) - E c (j +2 ) 



(g, J'l^intjc, j - l)(c, j - l|7€i nt |c, J + 2) | (g,j|Hint|g,j + 3)(g,j + 3|Hi nt |c,j + 2) 



(E e (j+2) - E gj)( E c(j-l) - E e (j+2)) 

2A b 



h 2 g 2 sA Q ^/(j + l)(j + 2) 



2A 2 h(2n + A b ) hfi(n + A b ) 



+ 



(E e (j+2) - Egj)(E g (j +3 ) - Egj) 
3oA b (2j + 3)(2A 2 + 9A b + 80 2 ) 
ft 2 2 (0 + A b ) 2 (20 + A 6 ) 



g 2 (2j + 3)V(j + l)(j + 2) a A o£ (Aft + 60) 2 
24frA 2 2 (A 2 + 50A b + 60 2 ) +v\a ) 



(g, jl^Intlg, 3 + l)(g,j + l|7frnt|g, 3 + 4 ) 

2 (E g (j+4) - Egj) 
(g,j\n Int \e,j + l)(e,j + l\H lnt \g, (j + 4)) 

(g,3\n Int \g,j + 3)(g, j + 3\H lnt \g, (j + 4)) 

2 (Eg(j+4) - E g j) 

(g,j|Hi nt |e, j + 3)(e, j + 3|7ft nt |g, j + 4) 



E g (j+i) ~ E sj 



E eU+i) ~ E sU+4) 



+ 



+ 



+ 



E e (j+1) ~ E gj 



E e(j+1) - E g(j+i) 



E sU+3) ~ E gj E g(j+3) ~ E g(j+4) . 



gVg + + W + 3jg + 4)aA 2 (A| - 30 2 ) 
8ftA 2 2 ( A 3 + f7A2 _ 9H 2 A 6 - 90 3 ) 



+ 0(a 2 



(g, j|Hint|g, j + l)(g, j + l\H lnt \e,j + 4) 

2(E (j +4 ) - -Egj) 

(g, j|7ft nt |e, j + l)(e, j + l\H Int \e,j + 4) 

2 (E e (j+4) - E gj) 
{g,j\n Int \g,j + 3)(g, j + 3\H lnt \e,j + 4) 

2 (E e (j+4) - E gj) 

(g,i|Wi nt |e,j + 3)(e,i + 3|Wi nt |e,i + 4) 



L^g(j+i) ^ Egj EgO+i) ^ E c(j - + 4) 



E c (j+1) - E gj E e(j+1) - E e (j+4) 



+ 



2 (E C 



0'+4) 



Egj) 



^g(j+3) _ ^g(j+3) ^ E e (j +4 ) 

1 1 



E c(j+3) - E gj E e(j+3) - E c(j+i) 



gVg + W + W + 3)(j + 4>A o£ (2A b + 5r>) 
6M7 2 (A* + 80A 3 , + 190 2 A 2 + 12ft 3 A 6 ) 



(e,j\Hi at \e,j + l)(e,j + l\H Iat \g, j + 4) 

2 (- E 'g(j+4) - Eej) 

(e,i|Wint|g, j + 3)<g, j + 3|Wi nt |g, j + 4) 



+ 



2{E i 



g(j+4) _ E e j) 



E c (j+1) - E cj E e(j+1) - E g(j+4) 



+ 



E g(j+3) - E cj E g(j+3) - E g(j+i) 



i|Wint|e, j + 3)(e, j + 3|Hi„ t |g, j + 4) (e,j|Wi nt |g,j + l)(g, j + l|Wi n t|g, 3 + 4) 



(Eg(j+4) - E e j)(E e 



+ 3) _ E e j) 



( E E(j+4) ~ E ej)( E g(3+l) ~ E g(3+4)) 



9 2 VU + W + 2)(j + 3)(j + 4>A e(5A fc - 120) 
12fiA 2 2 (A 2 - 70A b + 120 2 ) 



{e,j\n lnt \e,j + l)(e, j + l|Hi nt |e, j + 4) 

2 (E e (j+4) - E e j) 

<e,j|Hi nt |g,j + 3)(g, j + 3|?ft nt |c, j + 4) 

2 (E e (j+4) - E C j) 
(e, j\n Int \e, j + 3)(e, j + 3\H lnt \e }3 + 4) 

2 (E e (j+4) - E e j) 
(e, j I^Int |g, j + 1) (g, j + 1 |Hlnt [e, j + 4) 

2 (E e (j+4) - E j)(E g( j +1) - E e( j +4) ) 

gVg + W + 2 ) g + 3)(j + 4)«A 2 (A 2 - 30 2 ) 
8ftA 2 2 (Ag + £!A 2 - 90 2 A b - 90 3 ) 



+ 



E c(j+i) - E C j E e( j + i) - E c (j +4 ) 



Eg(j+3) - E 0J E g( j +3 ) - E c (j +4) _ 



+ 



E c (j +3 ) - E c j E e (j +3) - E e (j +4 ) 
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APPENDIX B: OSCILLATOR MATRIX ELEMENTS 



Here we give the explicit form of the functions L^o and L^o introduced in section llV Bl and derive the corresponding 
matrix elements y nm . The zeroth order contributions in the nonlincarity in section TlV Bl are denoted by: 



L 



LOl 



L L oi(g 2 ) 



-(g 2 



5 2 A 2 (2A b + 3ft) 



2ftA 2 (ft- 

4g 2 eA 



A h 



A 2 (A 2 + 3ftA b + 2ft 2 )^ 



L L oo+{g) = 
Lloi-(9 2 ) = 



5A0 



A b (ft + A b )' 

4g 2 eA 



A 2 ft(A b - 2ft)' 



The term independent of g is L^oi®) = — 3a/2M7. 
The terms linear in a and g are given by: 



L NO o+(a,g) 
LN02+(a 1 g) 
L N0 2-{a,g) 
LN02{a,g) 



3agA (A b + 2ft) 

~ hA b n(A b + ft) 2 ' 

3ag A (A 2 +6A b ft- 



4 Ml(A b 



3agA 



13ft 2 



ft) 2 (A 2 + 3A b ft)' 



?iA b (A b -3ft)(A b + ft)' 

4ae 



hA b tt 



Finally, the terms linear in a but quadratic in g are: 



LNOi s (j,a,g 2 
LNOi c (j,a,g z ) 
LNOi+{a,g 2 ) 
L NO i-{a,g 2 ) 
LN03(a,g 2 ) 
LN03+{a,g 2 ) 
L N 03-(a,g 2 ) 



6e 2 ag 2 3ag 2 A 2 ,[14(j + 1)A| - ft 2 A b (88 + 92j) - (3 + 5j)ftA 2 - (89j + 87)ft 3 ] 



~hA 2 w 



4?ift 2 A 2 (ft + A b ) 3 (A b - 3ft) 



2 _ 6e 2 ag 2 3«,g 2 A 2 [-14(j + 1)A 3 + (5j + 7)A 2 ft + ft 2 A b (92j + 96) + ft 3 (89j + 91)] 



~/iA 2 ft 3 
2ag 2 A e(4A 



29ftA 3 



4ftA 2 ft 2 (A b 
51ft 2 A? - 80A b ft 3 



■3ft)(A b 
- 124ft 4 ) 



-ft) E 



fift 2 (A b - 2ft)(A 3 + 3A 2 ft + 2A b ft 2 ) 2 
A 2 ft - 56ft 2 A b - 36A 3 ) 



6ag 2 eA (9A 3 1 A2 ° ccri2 



M7 2 (A 2 + 3ftA b + 2ft 2 )(A 2 - 2A b ft) ; 



a 5 2 A 2 (14A 



25A 2 ft 



130ft 2 A b - 261ft E 



8hA 2 n 2 (A b + ft) 2 (A 2 - 9ft 2 ) 
a 5 2 A e(A 3 + 3ftA 2 + 74ft 2 A b + 216ft 3 ) 



3^A 2 ft(A t 



-2ft) 2 (A 3 



a3 2 A e(24A 3 



■8A 2 ft + 19ft 2 A b + 12ft 3 )' 
i b - 936ft 3 ) 



3^A 2 ft 2 (A b - 2ft) 2 (A 2 - 7A b ft + 12ft 2 
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We can now give the expressions for y nm using rjj = 2 ^ J ^ , where < r]j < tt: 



V2j+l,2j+l 
V2j+l,2j+2 

U2j+l,2j+3 



2/2j+l,2j+4 — 



V2j+l,2j+5 



V2j+l,2j+e 



V2j+l,2j+7 



U2j+l,2j+S 



V2j+l,2j+9 
2/2j+l,2j+10 



-L N oo{j + l>a>ff) + L N oo{j,a,g) - cos^ [2L LO0 (.g) + L NO o(j + l,a,g) + L NO o{j,a,g)] (Bl) 
+ V.i + lsinr/j [L LO i-(g 2 ) + (j + l)L NO i- (a, g 2 )] , 
[2L LO o{g) + L NO a(j,a,g) + L NO o(j + 1, a, g)} sin^- 
+ Vj + Icosry^ [L i0 i-(5 2 ) + (j + l)L NO i- (a, g 2 )] , 

cos cos to V7T2 [l + (j + 2)iAr (a) + L LO i(g 2 ) + L N01g {j + 1, a, <? 2 )] , 
+ cos ^ sin to [iioo+(g) + L NO0 +(a, g){2(j + 1) + 1)] 
+ S i n !| cos to V(j + W + 2)L N02 -(<x, g) 

+ sin ^ sin to VJTT [1 + + l)Ljvo(a) - L LO i(g 2 ) + LjvoieC/, a, £ 2 )] , 
cos ^ cos to [iLoo+(.9) + iTO0+(«, ff)(2(j + 1) + 1)] 

- cos ^ sinte^j + 2 [1 + (j + 2)L NO (a) + L L01 {g 2 ) + L N01g (j + 1, a, <? 2 )] 



+ sin^-cos-^ [1 + (j + l)iivo(a) - L LO i(g 2 ) + L NO ic(j,a,g 2 )] 
- sin | sin to A /(j + l)(j + 2)L JV02 _( a , 5 ), 
cos ^ cos to y/(j + 2)(j + 3)L N0 2 (a, 5 ) 



v/(j + l)(j + 2)(j + 3)^03- (a, 5 2 ) 



+ cos -| sin VJT2 [i LO i+(.g 2 ) + (j + 2)LAr 01+ (a, g 2 

. 77, ?7 ? '+2 

+ sin — cos — — 

2 2 

_ sin ^ sin ^L NQ 2(a, < ? )v / (.7 + l)(j + 2), 

_ cos Hi si n to y/(j + 2)(j + 3)L N0 2 (a, <?) 

+ cos ^ cos to V7T2 [L LO i+(g 2 ) + (j + 2)L NO i+(a, g 2 )] 

_ sin | sin » A /(j + l)(j + 2)(j + 3)^03- (a, ff 2 ) 

_ sin | cos to iiV02 ( a7 ff ) v /(. ? + l)(j + 2), 

+ cos ^ costo v /( 7 + 2 )(. 7 +3)(. 7 +4) [i N03 (a, .g 2 ) - LAr (a)/2] 
+ cos Hi sin to V(j + 2)(j+3)Ljvo2+ (a, <7) 



+ sin | sin to v/(. 7 - + + 2)0' + 3) [-L W03 ( a , 5 2 ) - Ljvo(a)/2] , 

% • Vj+3 

— cos — sin — — 
2 2 

+ cos ^ cos to V(j + 2)(j + 3)ijvo2+ 5) 

+ sin ^ costo v /( j + i)( J+2 )(. 7 + 3) [-L N03 (a, g 2 ) - L NO {a)/2] , 



v /(. 7 + 2)(. 7 - + 3)(j+4) [W3(a,.g 2 ) - LAr (a)/2] 



cos | sin to v /(j + 2)(. 7 +3)(. 7 +4)L JV03+ (a, g 2 ), 

COS -f COS 

2 2 



y/(j + 2)(j + 3)(j + 4jL N0 3+(a, : g 2 ), 



and 
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V2j+2,2j+2 



2/2j+2,2j+3 



LNOo{j,&,g) - L NO o(j + 1,Q!,5) + ( 2L Loo(g) + L NO o(j,u,g) + L NO o(j + l,a,g))cosr)j (B2) 



2/2j+2,2j+4 



2/2i+2,2j+6 



V2j+2,2j+7 



V2j+2,2j+8 



V2j+2,2j+9 
V2j+2,2j+W 



SU1 cos-^v/j + 2 [1 + (j + 2)Ljvo (a) + L LO i{g 2 ) + L NO i g (j + l,a,g 2 )] 



- sin r?j + 1 [ l loi- (g 2 ) + (j + l)LNOi-(a,g 2 )] , 

• »?j »7i+i 

— cos — 

2 2 

- sin ^ sin 5^1 [L iO 0+(<?) + L NO0 +(a, <?)(2(j + 1) + 1)] 
+ cos ^ cos v /(. 7 + l)(. 7 + 2)i Ar02 _ (a, <?) 

+ cos ^ sin y/j + 1 [l + (j + l)L NO (a) - L LO i(g 2 ) + L NO ic(j,a,g 2 



sin -| sin V7+2 [l + (j + 2)L NO (a) + L LO i(g 2 ) + L NO i s (j + l,a,g 2 )] 



2/2j+2,2j+5 = 



— sin 


2 


cos 


Vj+i 
2 


— cos 


<h 


sin 


Vj+i 


2 


2 


+ cos 


<h 

2 


cos 


%+i 
2 


— sin 


'/.; 
2 


cos 


%+2 

2 


— sin 


r h 

2 


sin 


%+2 

2 


+ cos 


2 


cos 


2 


— cos 


»7j 


sin 


Vj+2 


2 


2 



sin | sin Hm^(j + 2 )(j + 3)L NQ2 (a, g) 



— sin 


2 


cos 


Vj+2 

2 


— cos 


<h 


sin 


V]+2 


2 


2 


— COS 


<h 

2 


cos 


Vj+2 

2 


— sin 


>h 

2 


cos 


2 


— sin 


Vj 
2 


sin 


2 


+ cos 


<h 

2 


sin 


%+3 

2 


+ sin 


'/.; 


sin 


Vj+3 


2 


2 


— sin 


2 


cos 


2 


+ cos 


% 
2 


cos 


»7i+3 

2 


— sin 


% 
2 


sin 


2 


— sin 


Vj 


cos 




2 


2 



v /(j + 2)(. 7 +3)0- + 4) [Ljvo3(a, g 2 ) - Ljvo(a)/2] 
V(j+2)(j+3)L 



v ^J+2XJ+3)(J+4)L A r 03+ (a, g 2 ) . 
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The matrix elements including the ground state are calculated separately because of its special form: 

2/oo = -2{L LO o(g) + L NO o{0, a,g)), (B3) 



Voi = cosy [l + L NO (a) +L LO i{g 2 ) + L NOlg (0,a,g 2 )] + sin y [L LO o+(g) + L NO a+{a,g)] , 

2/02 = _sin y C 1 + L No(a) + L LO i(g 2 ) + L NO i g (0, a,g 2 )] + cos y [L L oo+(g) + L NO o+{a, g)} 

2/03 = cos y %/2LAr 2(a, g) + sin y [L LO i+(g 2 ) + L NO i+{a, g 2 )] , 

2/04 = -siny\/2ijvo2(a,3) + cos — [L LO i+{g 2 ) + L W oi+(«, ff 2 )] , 

2/05 = cos y \/3 [Ljvo3 («, 5 2 ) - ^ivo (a) /2] + sin y V2Ltvo2+ (a, , 



2/06 
2/07 
2/08 



- sin y a/3 [iAr 03 (a, gt 2 ) - L NO {a)/2] + cos y \/2LAro2+(a, 9), 
sin Y^ L N03+(a, g 2 ), 
cos y V3L N0 3+ {a, g 2 ) ■ 



APPENDIX C: RATE COEFFICIENTS FOR THE OFF-DIAGONAL DENSITY MATRIX ELEMENTS 

We give the rate coefficients occurring in the Bloch-Redfield equation for the reduced density matrix, 

£01,01 = ^2/oo2/ii - ^2/ 2 o " ^yfi ~ 2^00,11, (CI) 

£02,02 = ^2/002/22 - ^2/o 2 o " ^2/ 2 2 2 - ^00,22, (C2) 

r ^ K K 2 K 2 ^ r ^ r fr«i\ 

M)3,03 — ^2/002/33 - ^2/00 ~ ^^33 ~ g 11,33 ~ 2 22,33 ' ^ ' 

^04,04 — -^2/002/44 - ^2/00 ~ ^^44 _ g 11,44 ~ 2 22,44 ' ^ ' 

£12,12 = ^2/112/22 - ^2/ 2 i - ^2/ 2 2 2 - 2-Ao.ti - ~£ 0,22, (C5) 
£13,13 = -^2/112/33 - ^pVii ~ ^g y 33 - 2^°° !l1 - 2^ 11,33 _ 2^ 22,33 ' 

£l4,14 = -^2/112/44 - Tgfll _ ^0^44 " T^ 00,11 - T^-WM ~ TJ^ 22 ' 44 ' 

r _ ^ K K 2 K 2 ^ r ^ r ^ r lr*Q\ 

^-23,23 — ^2/222/33 — ^p 22 ~~ ft/^ 33 ~~ 2 00,22 _ 2 11,33 _ 2 22,33 ' ^ ' 

£24,24 = ^2/222/44 ~ ^g2/22 - ^^44 _ T^ 00,22 _ 2^ 11,44 _ 2^ 22 ' 44 ' 

£34,34 = -^2/332/44 - ^2/33 ~~ ^^ 44 _ 2 11,33 _ 2 22,33 _ 2 11,44 _ 2 22,44 ' (CIO) 
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An,o2 = 7775(2/002/12 - 2/122/22) - ^(^02)^02^012/02 - G(w 12 )iVi2yiiyi2 
ftp 

-G(uJ 32 )N 32 y 13 y 23 - G(w 4 2)^422/l42/24, (Cll) 

£02,01 = 7^(2/002/12 - 2/122/11) - G(w 1)^0 12/0 12/02 - G{ui 2 i)N 21 y 22 y 12 

-G(w 3 i)iV3i2/i3y23 - G(w4l)^4l2/l42/24, (C12) 
^13. 23 = ^(2/332/12 - 2/I22/22) - G(wi 2 )iVl2(2/ll2/l2 - 2/122/33) 

ftp 

-G(to 02 )N 02 y m y m - G(uj 32 )N 32 y 13 y 23 - G(u> i2 ) N i2 2/142/24, (C13) 

£23,13 = 7^75(2/332/12 - 2/122/11) - G(w 2 i)iV 2 i(y222/i2 - 2/122/33) 
ftp 

-G(w i)iVoi2/oi2/02 - G(uj 3 i)N 31 y 13 y 23 - G(w 4 i)^4i2/i42/24, (C14) 

£14, 24 = 7^75(2/442/12 - 2/122/22) - G(wi 2 )iVi2(2/ii2/i2 - 2/122/44) 
ftp 

-G{uj 02 )N Q2 y 01 y 02 - G(uj 32 )N 32 y 13 y 23 - G{uj4 2 )N 42 yuy24, (C15) 

£24,14 = 7^75(2/442/12 - 2/122/11) - G(ui 2 i)N 2 i (2/222/12 - 2/122/44) 
ftp 

-G(w i)iVoi?;oi2/02 - G(w3i)-W3i2/i32/23 - G(w 4 i)7V 4 i?/i4?/2 4 . (CI 6) 

APPENDIX D: DIAGONAL REDUCED DENSITY MATRIX ELEMENTS 

The solutions of the FS A master equation (|53a| for the diagonal elements within the low temperature approximation 
(I5T|) reads: 

^OO(I) = CT + °U + °22 + °33 + a 



44 

r , , / £ll,33 £ll,44 

exp(-7r£ o,ii*) I °"n + °33 — 7. 7~r — 7. + °44 — 7. ~7 ~7 

i-00,11 + 1-11,33 + 1-22,33 — jS-00,11 + 1-11,44 + 1-22,44 



r -|0,0 £22,33 £22,44 

exp(-7r£ o,22£) cr 2 2 + 033 — 7 —7- —7 1" °44 



• exp(— 7r(£ii,33 + £22,33)^)^33 



£00,22 + £11,33 + £22,33 —£00,22 + £11,44 + £22,44 

£00,22 — £11,33 £11,33 



—£00,22 + £11,33 + £22,33 —£00,11 + £11,33 + £22,33 
£00,22 — £11,44 £11,44 



/ tr 1 r \A *-00,22 *-ll,44 . J-11,44 \ ,~ 1 

f exp(-7T(£n,44 + £22,44)t)ff 4 4 ~7 ±~7 ~7 + ~7 "7 ~7 ) , C°l) 

\ —'--00,22 + 1-11,44 + 1-22,44 —i-00,11 + 1-11,44 + 1-22,44 , 

en (t) = -cxp(-7r£ 00ill i)cr5 ) 1 



.0 




£ll,33 




33 


— £00,11 


+ £ll,33 " 


1" £22,33 


.0 




£ll,44 




44 


— £00,11 


+ £ll,44 " 


\~ £22,44 



- exp(-7r(£ o,ii + £n,44 + C 22A i)t)al 4 — — , (D2) 

ff22(0 = - exp(-7r£ o,22i)o-22 



exp(-7r(£ o,22 + £11,33 + £22,33)i)o~3; 



33 " 



£22,33 

-£00,22 + £11,33 + £22,33 



£22 44 

exp(-7r(£ o,22 + £11,44 + £22,44)^)0-44—7 —7^ —7 > ( D3 ) 

— 1-00,22 + 1-11,44 + 1-22,44 



24 

o-aa (t) = exp(-7r(£n,33 + £22,33)^33' ( D4 ) 

0-44 (*) = exp(-7r(£n,44 + £22,44)t) cr 44- (D5) 

I 
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